function z = mult (x, y, d1, d2)

   ## usage:  z = mult (x, y, d1, d2)
   ##
   ## nd-array multipication, along dims d1, d2

   dx = size(x) ; dy = size(y) ;

   if (nargin < 3)
      d1 = length(dx) ;
      d2 = 1 ;
   endif

   I1 = [dx(1:d1-1) dx(d1+1:end)] ;
   I2 = [dy(1:d2-1) dy(d2+1:end)] ;

   xx = reshape(x, prod(I1), dx(d1)) ;
   yy = reshape(y, dy(d2), prod(I2)) ;

   zz = nanmult(xx, yy) ;
   # zz = xx * yy ;
   z = reshape(zz, [I1 I2]) ;

endfunction
